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Distributed and Recursive Parameter Estimation 
in Parametrized Linear State-Space Models 

S. Sundhar Ram, V. V. Veeravalli, and A. Nedic 
Abstract 

We consider a network of sensors deployed to sense a spatio-temporal field and estimate a parameter 
of interest. We are interested in the case where the temporal process sensed by each sensor can be modeled 
as a state-space process that is perturbed by random noise and parametrized by an unknown parameter. 
To estimate the unknown parameter from the measurements that the sensors sequentially collect, we 
propose a distributed and recursive estimation algorithm, which we refer to as the incremental recursive 
prediction error algorithm. This algorithm has the distributed property of incremental gradient algorithms 
and the on-line property of recursive prediction error algorithms. We study the convergence behavior of 
the algorithm and provide sufficient conditions for its convergence. Our convergence result is rather 
general and contains as special cases the known convergence results for the incremental versions of the 
least-mean square algorithm. Finally, we use the algorithm developed in this paper to identify the source 
of a gas-leak (diffusing source) in a closed warehouse and also report some numerical results. 

I. Introduction 

A sensor network consists of sensors that are spatially deployed to make observations about a process or 
field of interest. If the process has a temporal variation, the sensors also obtain observations sequentially in 
time. An important problem in such networks is to use the spatially and temporally diverse measurements 
collected by the sensors locally to estimate something of interest about the process. This estimation 
activity could either be the network's main objective, or could be an intermediate step such as in control 
applications where the sensors are also coupled with actuators. 
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third author is with the Dept. of Industrial and Enterprise Systems Engg., University of Illinois at Urbana-Champaign. They 
can be contacted at {ssrinivS, vvv, angelia}@uiuc . edu. This work was partially supported by a research grant from 
Vodafone and partially by a research grant from Motorola. 
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In this paper, we consider a parameter estimation problem when each individual sensor observation 
process can be modeled as a linear state-space process that is parametrized by an unknown parameter 
of interest, and also perturbed by process and observation noise. State-space models arise directly, or as 
linear approximations to non-linear models, in many applications. As an example, we will later discuss 
the problem of localizing a gas-leak in a warehouse. 

We propose a distributed and recursive estimation procedure, which is suitable for in-network pro- 
cessing. Each sensor locally processes its own data and shares only a summary in each time slot. The 
sensors form a cycle and update incrementally, whereby each sensor updates the estimate using its local 
information and the received estimate from its upstream neighbor, and passes the updated estimate to its 
downstream neighbor. In this way, there is a reduction in communication within the network at the cost of 
increased local sensor processing. This can significantly reduce the total network energy used, especially, 
when the sensors communicate over a wireless medium. Furthermore, the sensor updates are generated 
recursively from every new measurement using only a summary statistic of the past measurements. This 
has two benefits. Firstly, the network has a (possibly coarse) estimate at all times, which is important in 
applications that require the network to react immediately to a stimulus and make decisions on-line. For 
example, in a network that is deployed to monitor gas leaks the network should raise an alert depending 
on the level of the leak intensity. An additional benefit is that each sensor can purge its old measurements 
periodically since only a summary of constant size is used to update the estimates. This can significantly 
reduce the memory requirements at the sensors. 

Our approach is in contrast with traditional estimation methods such as the maximum likelihood and 
least-squares which are centralized, i.e., the measurements collected by the spatially distributed sensors 
are routed through the network to a single location (fusion center) where estimates are computed. In this 
case, the network energy is mainly consumed in routing the measurements to the fusion center, which 
can be inefficient in terms of energy consumption. The problem of centralized recursive estimation in 
linear state-space models is an old problem in system identification. We refer the interested reader to [1] 
for a survey of these methods for linear state-space models. The problem has also generated considerable 
interest in the neural networks community where the EM algorithm is used as a tool to learn the parameters 
[2]. A related algorithm is the parallel recursive prediction error algorithm proposed in [3] that updates 
the components of the parameter vector in parallel. 

The literature on distributed estimation is somewhat limited. A distributed maximum-likelihood algo- 
rithm is discussed in [4] and a distributed expectation-maximization algorithm is discussed in [5]. In 
[6], the incremental (sub)gradient algorithms of [7] are used to obtain distributed least-square estimators. 
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Distributed linear least-squares are discussed in [8] without an explicit point-to-point message routing. 
All of these algorithms are distributed but not recursive. 

This paper extends our earlier work [9], where we considered the problem of recursive and distributed 
estimation for stationary models. To the best of our knowledge, there is only one other related study [10] 
that deals with both distributed and recursive estimation. There, incremental versions of the least-mean 
square algorithm and the recursive least-squares are developed to solve the linear least-squares problem. 
In both studies [9] and [10], the models are not auto-regressive. 

Our contribution in this paper is the development and convergence analysis of a general distributed 
recursive algorithm for parameter estimation in parametrized state-space models. Our results are more 
general than those of [10], which follow as a special case. 

The rest of the paper is organized as follows. We formulate the problem in Section|IIl and then introduce 
our notation in Section |llll We give an overview of the algorithm in Section |IVl We then discuss the 
standard recursive prediction error algorithm algorithm [1] and the incremental gradient algorithm of [7] 
in Section lYl These algorithms are at the heart of our distributed algorithm presented in Section |Vll where 
we also state our main convergence result. We prove the convergence of the algorithm in Appendix [A] 
We discuss some simple extensions in Section IVlIl We report some experimental results obtained by our 
method as employed to localize the sourse in a gas leak problem in Section IVIIII We conclude in Section 
US 

II. Problem Formulation 

We consider a network of m sensors, indexed 1, . . . ,m, deployed to sense a spatio-temporal diverse 
field to determine the value of some quantity of interest, denoted by x, with x G Jft''. We denote the true 
value of the parameter by x*. 

We assume that time is slotted and each sensor sequentially senses the field once in every time slot. 
We model the measurement sequence of sensor z as a random process {Ri{k;x)} with the following 
dynamics 

Q,{k + l;x) = D,{x)Qi{k;x) + Wi{k;x), 

Ri{k + 1; x) = HiOiik + l;x) + Vi{k + 1). (1) 

Here, {Wi{k;x)} is the process noise, {Vi{k)} is the measurement noise. Hi is the observation matrix 
and Vi{k + 1) is the measurement noise of sensor i. The process {@i{k;x)} can be interpreted as the 
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temporal process obtained by sampling a spatio-temporal diverse field at the location of sensor i. At this 
point, we do not assume any knowledge on the joint statistics of Qi{k + 1; x) and Qj{k + 1; x)Q 

We denote by ri{k) the actual measurement collected by sensor i at time slot k, i.e., ri{k) is a 
realization of Ri{k;x*). The processes {Wi{k;x)} and {Vi{k)} are zero-mean i.i.d. random sequences. 
The quantities Di{x), Hi, Cov(Wj(/c; x)) and Coy{Vi{k)) are available only at sensor i. Moreover, at all 
sensors a set X is available that satisfies the following properties 

1) The set X is closed and convex; 

2) The true parameter x* is contained in the set X; 

3) The system in ([T]) is stable, observable and controllable for all x ^ X. 

Note that X may even be the entire SR'^. The problem is to estimate the parameter x from the collection 
of sensor measurements {ri{k)} with an algorithm that is: 

1) Distributed: Sensor i does not share its raw measurements {ri{k)} with any other sensor. 

2) Recursive: At all times, sensor i stores only a summary statistic of a constant size, i.e., the size of 
the statistic does not increase with the number of measurements collected by the sensor. 

III. Notation 

All the random variables are defined on the same probability space T = {Q,!F,V). If a; G 17, is the 
outcome of an experiment, then for a random process {Y{k;x)} that is parametrized by x, we define 
y{k) =: Y^{k; x*), i.e., y{k) is the value of the random variable Y{k; x*) corresponding to the outcome lo. 
According to this notation, ri{k) and 9i{k) are the realizations of Ri{k; x*) and @i{k; x*) that correspond 
to the same outcome lo. 

We let I denote the set of sensors, i.e., T := {1, . . . , m}. Further, we assume that &i{k; x) and Ri{k; x) 
are vectors of dimensions q and p, respectively. They are the same for all sensors i G xj^ We write 
R^{x) to denote the collection of random variables x) , . . . , Ri{k; x)}, which should be viewed as 

a collection of random variables parametrized by x and not as a function of x. Furthermore, in line with 
our notation, rf denotes the realization of R^{x*), i.e., rf denotes the collection {rj(l), . . . ,ri{k)}. 



'Even if some information is available we ignore it. This aspect is discussed in detail in Section IVII-BI 
^We make this assumption only for the sake of clarity. Our analysis applies to the general case where the dimensions q and 
p can be sensor dependent. 
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IV. Algorithm Overview 

A standard estimation procedure defines the estimate as the minimum of a suitably defined cost that is a 
function of the observations and the unknown parameter. For example, the maximum likelihood estimator 
minimizes the negative of the log-likelihood function. The form of the cost function determines whether 
there is a distributed and recursive minimization procedure. Further, the cost function also determines 
other properties of the estimator such as unbiasedness, consistency, minimum variance etc. (see [11]). 

Except in some very special estimation problems, it is impossible to find a cost function that supports 
even a centralized recursive procedure and also generates a 'good' estimate. In this paper, we develop a 
distributed and recursive estimator that is only consistent, i.e., the estimate converges to the correct value 
X* as the number of available measurements becomes infinitqj. The estimates are biased but this is the 
price that is to be paid to obtain a distributed and recursive procedure. Thus, there are two aspects to the 
problem. The first is to choose a suitable cost function, and the second is to develop a distributed and 
recursive minimization procedure. We will first discuss the cost function that is used, and then give an 
overview of the minimization algorithm. 

A. Cost function 

Suppose that each sensor has made N measurements and we want to estimate the parameter x from 
these measurements. As mentioned, the cost is a function of both the available measurements and the 
unknown parameter. Therefore, we denote the cost function as fN{x;r^). 

For X G X, we assumed that the system in ([T|) is stable, observable and controllable. The Kalman gain 
for the system therefore converges to a finite time-invariant value [12]. Let Gi{x) be the Kalman gain 
for the state-space system in ([T]), which is determined from Di{x),Hi, Cov {Wi{k] x)) , and Coy{Vi{k)) 
as the solution to the Riccati equation [1]. Using Gi{x) define 

(/.i,fc+i(x;rf) = {Di{x) - Gi{x)) (l)i,k{x;r^-^) 
+ Gi{x)ri{k), 

9i,k+i{x;ri) =i?i(/)i_fc+i(x;rf), (2) 

with r^') = fii{x). Observe that gi^k+i{x',r^) is linear in rf for each x. Furthermore, for any 

X E 3f?'^, gi,k+i{x', r^) viewed as a function of rf is an one-step prediction filter (henceforth, referred to 



'This statement is technically imprecise and will be clarified later. 
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as a predictor) for the random process {Ri{k + l;x*)}. Thus, {gi^k+i{x\r^)}x^-)i^d is a predictor family 
parametrized by x. 

We will choose our cost function to be 

N m 



^ 7V rn 



k=l i=l 



(3) 



and our estimator to be 



xn = arg min /7v(a;; r^). 

We next provide an intuitive explanation as to why this choice of cost function should generate consistent 
estimates. First, note that the vector gi^k+i{x*;r^) is the steady-state Kalman predictor for Ri{k + 1; x*) 
since {ri{k)} is a sample path of the random process {Ri{k; x*)}. Among other properties, the steady 
state Kalman filter is asymptotically optimal in a mean square sense in the class of linear time-invariant 
predictors. Thus, x* will minimize 



N m 



f(x) = lim — We 

k=l i=l 

= E[Mx;R^{x*))] 



R^{k;x*)-gi,k{x;RtHx*)) 



Here, the expectations are taken with respect to the random process {R^{x*)}. Therefore, one can expect 
that xj\f, which is the minimum of fj\f{x;r^), might in the limit be equal to x*, since x* is the minimum 
of fix), the Umit of E[fN {x;R^{x*))] . 

B. Overview of the minimization procedure 
Observe that fN{x,r^) can be written as 

m 

where 



i=l 



1 ^ II 

fi,N {x; rf) =-j^Yl Ir ^(^) " 3i,ki^-^ ^) 



k=l 

Suppose we are interested in a non-recursive but distributed solution the problem. Then, the incremental 
gradient algorithm [7] can be used to minimize fj\i{x;r^). In each time slot, the algorithm cycles the 
estimate through the sensor network. Sensor i receives the estimate Zi^i^k+i from sensor i — 1 at time 
slot A; + 1, and generates a new estimate Zj.fc+i using V/j(zj_i^fc+i). The new estimate is then passed to 
sensor i + 1 for i < m, and to sensor 1 for i = ?n and, thus, the estimate is cycled through the network 
for each sensor to update. An illustration is shown in Fig.[T] 
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Fig. 1. A network of 10 sensors with incremental processing. The estimate is cycled through the network. The quantity Zi^^+i 
is the intermediate value after the sensor i update at time k + 1. 



Now consider the complementary problem, i.e., a centralized but recursive solution. By recursive 
we mean, the algorithm should be able to obtain xn+i directly from xn, the new measurements 
ri 7V+1; • • • , '^m.TV+i, and some summary statistic of the past observations. This is generally not possible. 
Nevertheless, it is possible to use the recursive prediction error algorithm of [1] to obtain recursive 
approximations to {xk} such that the approximate sequence converges to the same limit as {xk}. Thus, 
while we do not minimize the chosen cost function at each step, the new sequence {xk} is still 
consistent. 

The algorithm that we propose is a combination of the incremental gradient algorithm and the recursive 
prediction error algorithm. We therefore refer to the algorithm developed in this paper as the incremental 
recursive prediction error algorithm. 

V. Preliminaries 

We first evaluate some quantities that will be useful to us in the analysis. To make the paper self- 
contained, we then discuss the incremental gradient method [7] and the recursive prediction error algorithm 
[1] in this section. 

A. Some more notation 

For later reference, we obtain the form of the gradient of the predictor Qi^k+i {x', rf) . Define Fi{x) = 
Di{x) — Gi{x)Hi and for convenience rewrite Q as follows 

(/>i,fc+i(2;;rf) =Fi{x)4)i^k{x]r^i~^) + Gi{x)ri{k), 

5i,fc+i(2;;rf) =Hi(j)i^k+i{x;r^). (4) 
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Let x^^) denote the £-th component of x, and define 



'k,k 



dF.jx) 
dG^jx) 

dxW ■ 



Thus the gradient Vgi^kix'^r'l ^) is the p x d matrix, 



(5) 



By differentiating in we can immediately see that 



c2+i(-;rf) 

9Lk+i{x;r'y) 





vWf,(x) F,(x) 





H, 



(PLk+i{x;r^ J 



k\ 



+ 



Gi{x) 
VWG,(x) 



riik), 



(6) 



B. Incremental gradient descent algorithm 

For differentiable optimization problem of the form 



mm 

x6X 



i=l 



the standard gradient descent method, with projections, generates iterates according to the following rule: 



Xk+l 



Vx 



Xk - Ofc+l 



Ev/,(xfc) 



Here, the scalar Uk+i > is the step-size and Vx denotes the projection onto the set X. This method 
is centralized in the sense that it requires the gradient information of each fi{x) at the current iterate Xk 
in order to generate the new iterate Xk+i- In our setting, however, the gradient information V fi{x) is 
distributed since fi{x) is known only locally at sensor i. Thus, the standard gradient descent method is 
not adequate. 

To deal with the distributed nature of the sensor network information, we consider the incremental 
gradient method to minimize f{x), without the sensors explicitly sharing the functions fi{x) (see, [7], 
[13] and the references therein). In this algorithm, the iterates are generated according to 



Xk — ^m,k — ^0,fc+l) 
Zi,k+1 = Vx [Zi-l,k - Ok+l^ Mzi^i^k)] 



(V) 



The key difference between the standard gradient and incremental gradient method is that the standard 
gradient method generates iterates by using the gradient information of all functions fi{x) at the same 
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(current) estimate Xk, while the incremental method generates iterates through a cycle of intermittent 
adjustments Zi_i^k+i using only one function at a time, i.e., the gradient V/j(zj_i fc+i), so that all 
functions /j are processed within a cycle (see Fig.[T]for an illustration). The convergence of the incremental 
gradient method has been studied in [13], [7], [14] under different assumptions on the functions fi{x) 
and the step-size rules. 



C. Recursive prediction error algorithm 

Here, we discuss the standard recursive prediction error algorithm (RPE) for a parameter estimation 
problem (see [1]). To avoid confusion with the notation in the rest of the paper, we suppress the subscript 
i and consider the problem of estimating x from observations of a random process {R{k;x)} with the 
following dynamics: 



@{k + 1; x) = D{x)Q{k; x) + W{k; x) 
R{k + l;x) = He{k + l;x) + V{k). 



(8) 



The RPE algorithm generates estimates of x by applying suitable approximations to the iterates generated 
by the gradient descent method as employed to solve an appropriate optimization problem. In particular, 
on the set X, the true parameter value x* minimizes the following function: 

2" 



lim 



1 ^ 



k=l 



(9) 



\R{k;x*)-9k{x;r''-'{x*)) 
When the standard gradient descent method is used to minimize /, the iterates x^+i are given by 

Xk+i = Vx [xk - afc+iV/(xfc)] . (10) 

The RPE algorithm obtains a sequence {x^} by using two approximations to the sequence {x^}. 

The gradient of f{x) is not available in ([TOb . but instead the sequence {r{k)} is available. The first 
approximation is a least mean-square (LMS) type approximation replacing the actual gradient V/(xfc) 
with an empirical gradient. Let us denote the iterate sequence corresponding to this approximation by 
{xfc}, for which we have 



Xk+i = Vx 



fe+l^ 



where 



V/fc+i(xfc;r 



Xk - Cik+i\7 fk+i{xk;r 



2 [Vgk+i{xk\r^)] {r{k + l)- gk+i{xk\r^)). 
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The gradient V(7fc+i(x; r'^), can be obtained from (O and the extended representation of 7]f^-^^{x;r^) in 
The problem is that even with this approximation the sequence {x^} cannot be obtained recursively. 
Observe that to exactly evaluate gk+i{xk',r^) and Vgk+i{xk;r^) one would need the entire vector r^. 
To accommodate the recursive computations, we use another approximation 



(/)fc+i(xfc;r'^') 



F{xk) 
vWF(xfe) F{xk) 





G{xk) 


+ 


_ vWG(xfc) _ 



r{k), 



H 
H 



(11) 



(/)fc+i(xfc;r'') 

cSi(^fc;^'') 

Changing notation to reflect the approximations and re-ordering the equations, the resulting RPE algorithm 
can be written as follows, for £ = 1, . . . , d, 



hk+i 




H 







V'fc+i 


?/c+l 







H 







efc+1 

^k+l 
^k+l 
Xk+1 

Xk+2 



r{k + l) -h 



k+l, 



■ ^k+1 



Cfc+l, 



[Xk+l] , 

vWF(xfc+i) F{xk+i) 



i^k+i 
(^) 



+ 



G{xk+i) 
vWG(xfc+i) 



r(A; + l). (12) 



The algorithm is initialized with values for ipi,Xi ^^'^ ^o- Observe that to update Xk the algorithm 
requires only r{k + 1), xi+i) • • • ) xi+i ^^d ipk+i, and therefore, it is recursive. 

In summary, the iterates of the RPE algorithm are obtained from the standard gradient descent iterates 
with the following two approximations: 

1) An LMS-like approximation for the gradient, and 

2) An approximation to make the LMS approximations recursive. 

The following theorem provides some sufficient conditions guaranteeing that the iterates generated by 
the RPE algorithm asymptotically converge to a minimum of f{x). The theorem is based on the results 
from [1]. 

Theorem 1: Let the following conditions hold. 

1) The set X is a closed and convex set containing x* . Furthermore, the system in (IS]) is stable, 
observable and controllable for all x ^ X. 
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2) The matrices F{x) and G{x) are twice dijferentiable for all x & X. 

3) The fourth moments ofV{k) are bounded. The second moments ofW{k;x*) are bounded. 
Moreover, let the step-size ak be such that kak for some positive scalar ^. Then, the iterates Xk 
generated by the RPE in di2D converge to a local minimum of f{x) in ((91) over the set X, with probability 
1. 

Theorem [T] follows from Theorem 4.3 on page 182 and the discussions in pages 172 and 184 of [1]. 
The conditions for convergence of the algorithm are extremely weak. Note that the algorithm guarantees 
convergence only to a local minima and not necessarily to the global minimum x* of f{x). Of course, 
when the function f{x) is convex this implies convergence to a global minimum. 

VI. Incremental Recursive Prediction Error Algorithm 
As discussed in Section |IV-A[ when there are multiple sensors, the true parameter x* minimizes 



f{x) = y lim E 

i=l k=l 
m 



\Ri 



(13) 



i=l 



We combine the incremental gradient algorithm in ^ with the RPE algorithm in (fT2l ) to develop an 
incremental recursive prediction error (IRPE) algorithm. The main idea of the IRPE is to use an RPE 
like approximation for the gradient term in the incremental gradient algorithm d?]). Formally, the iterates 
are generated according to the following relations for i G T, and I = 1, . . . , d, 



Xk 



hi^k+i 




Hi 







■ipi,k+i 


?i,fc+l 







H. _ 




M,k+1 



ej,fc+i 

^i,k+l 



^i,k+l 



Zi,k+1 



ri{k + 1) - hi^k+1, 



'l,fc+l - "fc+l [^i,k+l ) ^hk+l, 
T 



(1) 



Sd) 
^i,k+l 



i'i,k+2 








X.i,k+2 





fc+1 
'Px\zi,k+i\, 

Fi{zi^k+i) 

vWFi(z,,fc+i) Fi{z,,k+i) 



Xi^k+1 



+ 



Gi{Zi,k+l) 

V^'^Gi{zi,k+i] 



(14) 

(15) 
(16) 
(17) 
(18) 

ri{k + l). (19) 



The initial values for the recursion are fixed at xo = Xg, V'j.i = ipi,s and xfl = xfl- To see that the 



algorithm has a distributed and recursive implementation assume sensor i — 1 communicates Zi^i^k+i to 
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sensor i in slot k + 1. Sensor i then use^ ri{k + 1) to updates the iterate zi^i^k+i to generate Zi^k+i- 
This is then passed to the next sensor in the cycle. Observe that in updating fc+i, sensor i requires 
only x^ik+v ■ ■ ■ XilVi V'i.fc+ii which were calculated by sensor i in the previous time slot. Thus, the 
algorithm is recursive and distributed. Furthermore, note that sensor i only needs to know its own system 
matrices Hi,Fi{x) and Gi{x). 



A. Convergence result 

The iterates generated by the IRPE method are three approximations away from the iterates generated 
by the standard gradient descent method. The first approximation is in going from the standard gradient 
algorithm to the incremental gradient algorithm, the second is in approximating the gradient of the function 
with an LMS-like empirical gradient and the third is in calculating the empirical gradient recursively. 
Therefore, it is not clear if the iterates will converge to x* . We next state a theorem that provides sufficient 
conditions for the convergence of the iterates generated by the IRPE algorithm. 

Theorem 2: For all i ^Z, let the following conditions hold 

1) The set X is a closed and convex set containing x* . Furthermore, the system in (jSj) is stable, 
observable and controllable for all x £ X. 

2) The matrices Fi{x) and Gi{x) are twice differentiable for all x £ X. 

3) The fourth moments ofVi{k) are bounded. The second moments of Wi{k; x*) are bounded. 
Moreover, let the step-size be such that ka^ /i for some positive scalar ^. Then, the iterates x^ 
generated by the IRPE algorithm in di4D-fl79l) converge to a local minimum of f{x) in l \13i over the set 
X, with probability 1. 

Note that the result implies that for each i the iterates Zi^k+i converge to the same local minimum. 
Thus the algorithm is not necessarily consistent. 

There is alternative way to interpret the IRPE algorithm. In particular, consider a centralized scheme 
where the sensors immediately communicate their measurements to a fusion center. Now, at the fusion 
center, the RPE algorithm can be used to estimate the parameter x. For the specific model of our interest, 
there is a hidden structure in the RPE algorithm that permits an incremental implementation. Thus, the 
IRPE algorithm can also be viewed as an incremental implementation of a centralized RPE algorithm. 
Since this hidden structure in the RPE algorithm is not easily identified, we have not used this alternative 



We are assuming that sensor i obtains its measurement before it receives the iterate. From an implementation perspective, 
each time slot can be divided into two parts. In the first part, the sensors make measurements and in the second part they process. 
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approach to actually present the algorithm. However, we use this approach to prove the convergence 
result stated in Theorem [2l The proof is provided in Appendix [A] 

B. Communication requirements 

Incremental algorithm can potentially require less communication than centralized schemes. In a 
centralized scheme, in every slot each sensor has to communicate its measurements to a fusion center that 
is 0(1) meters away on average. Summed over the m sensors in the network, the total communication 
requirement in a centralized scheme is 0{m) bit meters per slot. In the incremental scheme, each sensor 



[6]. Therefore, the total communication required in the incremental processing scheme is O(y^logm) 
bit meters per slot. 

C. Centralized versus incremental: Tradeoff 

In our analysis we do not use any information about the joint statistics of the random process {Qi{k; x)} 
and {Qj{k;x)}. When this is the case the performance of the IRPE is identical to the performance of 
the centralized RPE algorithm. 

Suppose some information about the join statistics is available. This information cannot be used in a 
distributed system because at most one sensor's measurement is known at a single location at any time. 
Thus, the joint distribution information is not useful to the RPE. 

A centralized system, on the other hand, can potentially use the joint density information to obtain a 
cost function f{x) that generates estimates with better properties. As an example, suppose that Qi{k; x) = 
Qj{k;x) for all /c > 1 and i,j G X, which corresponds to the case when all sensors sense a field with 
no spatial variation, synchronously at time mk. Define H, respectively V{k + 1), to be the block matrix 
obtained by stacking the matrices Hi, . . . , Hm, respectively vectors Vi{k + 1), . . . , Vm{k + 1). Then, the 
centralized measurements R{k; x) have the following evolution: 




meters away on average, as discussed in 



G(A; + l;x) = D{x)Q{k- x) + W{k; x) 



R{k + l]x) = HQ{k + l;x) + V{k + 1). 



(20) 



The corresponding time-invariant Kahnan predictor is given by 



</>fc+i(x; r'^) = {D{x) - G{x)H) r 



) + G{x)r{k), 



9k+i{x]r 
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Notice that the predictor gi^k+i{x; r^) for the {k + l)-st measurement of sensor i is a function of the past 
measurements made by sensor j, j i. Using this predictor we can define a cost function in a manner 
similar to ([3]l. As we will see in Section [Villi the nature of the cost function may be significantly better 
in terms of the number of local minima and the RPE applied to the system in (l20l) may have a better 
performance. 

To summarize, there is an implicit tradeoff when we use the IRPE. Potentially better estimates may 
be obtained by a centralized scheme when the joint statistics of the process 0j(A;;x) and Qj{k\x) are 
available. This is indicated by our numerical results in Section IVIIII 

VII. Extensions 
We next discuss some extensions to the IRPE algorithm. 

A. Hybrid scheme 

Let us consider an alternative network architecture where the network of m sensors, divided into vtic 
clusters of approximately equal size, is deployed in a unit square. Each cluster has a cluster head that 
is a neighbor to all the sensors in the cluster. We can develop a hybrid algorithm that is centralized 
intra-cluster and distributed inter-cluster. Each cluster head collects all the measurements made by the 
sensors in the cluster, and then the cluster heads use the IRPE algorithm to estimate x without sharing 
their measurements. 

Note that, as each sensor is in the neighborhood of its cluster head, it is still required to only 
communicate to a neighbor. The cluster heads might have to communicate over larger distances. The 
total inter-cluster communication is 0{mc) bits per meter and the total communication in a cluster 

bits over an average distance of Therefore, the total communication is 
^ (y'^^^°^ (™")) ^^^^ meter. The benefit is that the cluster heads can use any information 

that is available about the joint statistics of the processes seen by the sensors in the cluster. 

B. Distributed and recursive regression 

In the problem we have studied, the actual sensor measurement sequence {ri{k)} is a sample path of 
the random process {Ri{k; x)} for x = x* . While we did not assume to know the value of x*, we did 
assume that for some x ^ X the actual system is correctly modeled by ([T|). 

In practice, it is very difficult to obtain the correct description and often approximate models are used. 
In the context of our problem, this means that ([1]) need not necessarily be the correct description of the 
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actual system dynamics for any value of a: G X. The minimum of f{x) is now interpreted as the value 
of X that chooses the state-space system that best approximates the actual system among all the systems 
generated as x ranges over X. 

Theorem 4.3 of [1] concludes that even in this case, under some weak regularity conditions on the 
actual measurement sequence (instead of Condition 3) the iterates generated by the RPE still converge 
to a local minimum of f{x). Since the IRPE was proved to be equivalent to the centralized RPE, the 
above statement extends to the IRPE algorithm also. 

C. Other extensions 

We have not included an explicit input in modeling the system. Much of the analysis immediately 
follows where there is a deterministic open-loop input {ui{k)} that drives the system in ([Til. Of course, 
{ui{k)} should be known to sensor i. Another immediate extension is to the case when the matrix Hi 
and noise Vi{k) are also be parametrized by x. Finally, we remark that rate of convergence results are 
available for the RPE through a central limit theorem and these can be extended to the IRPE. 



We next consider a gas-leak problem to illustrate the concepts developed in the paper. We assume 
that a wireless sensor network is deployed inside a warehouse where gas tanks are stored. The network 
objective is to localize a leak, when one occurs. We use a two-dimensional model for this scenario, which 
is appropriate when the gas is significantly heavier than air. In any case, the extension to three dimension 
is immediate. We also remark that we have used the gas leak problem as only a representative example; 
the analysis is more generally applicable to heat and other diffusing sources. 

Leak model: We asume that the leak occurs at time t = and that the network detects the leak 
immediately. We model the leak as a point source at x = {xi, X2). Each sensor sampling has a duration 
of 1 time unit. The leak intensity is modeled as a piece-wise constant function, i.e., the leak intensity is 
equal to Ik during the time interval [A; — 1, fc) for A; > 1. Across sampling intervals, the leak intensity 
values vary according to the following Markov process: 



Here, p is a known scalar and {^(A:)} is a sequence of i.i.d. Gaussian random variables with zero mean 
and variance cr^. Thus, the intensity evolves in time as follows: 



VIII. Application 



Iik + 1) =pI{k) + S{k). 



(21) 



00 




(22) 



k=0 
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where rect(t) is the rectangular function taking value 1 in the interval [0, 1] and zero elsewhere. 

Medium model: We model the warehouse as a rectangular region with known dimensions li x I2. 
Without loss of generality, we let the warehouse to be the region D = [0, li] x [0, 12], and we denote the 
boundary of the warehouse by dD. 

The medium is characterized by the diffusion coefficient of the gas, boundary conditions and initial 
conditions. We use C{y, t; x) to denote the concentration at a point y at time t when the source is at x. 
We make the following assumptions. 

1) The diffusion coefficient of the gas is the same everywhere in the warehouse. We use v to denote 
this value. 

2) The boundaries of the room are insluated, i.e., there is no leakage out of the room, i.e., ^"^g^ '^^ = 
0,Vs G dD. 

3) At time t = the concentration is everywhere in the room, i.e., C(-,0;x) = 0. 
Observation model: Let Sj be the location of the i-th sensor. We assume that all sensors sense at the 

beginning of each time slot, i.e., at time k. Then 

R^{k;x)=C{si,k;x)+Ni{k), (23) 

where Ni{k) is a zero mean i.i.d. measurement noise with known variance cr^. 

Transport model: We assume that the transport of the gas in the warehouse obeys the diffusion equation. 
Therefore, 

= uV^C{y,t;x)+I{t)6{y-x), (24) 



dt 

with the initial and boundary conditions 

C{s, 0; x) =0 for all s e D, 

=0 for alH > and s e 

dt 

Here, is the Laplacian differential operator and S is the Dirac delta function. 

A. Problem statement and related literature 

The medium characteristics are completely known, i.e., li, I2, and u are known. The sensors' sampling 
duration and the measurement noise variance cr^ are also known. Further, the variance tr^ of S{k) is 
known. The problem is to determine the location of the point source x from the sensor measurements in 
a distributed and recursive manner. To solve the above problem we first show that, as a consequence of 
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the assumptions that have been made, the sensor measurements follow a state-space model. We then use 
the IRPE algorithm developed in the previous sections to solve the problem. 

We next compare and contrast the models described above with the models used in literature. The point 
source model is a common model for diffusing sources and has been extensively used in localization 
studies [15]-[18]. The random time-varying source intensity model is more realistic compared to the 
constant intensity [15], [17], [18] and instantaneous intensity models that are usually studied. Localization 
of sources with time-varying intensity have been studied in a centralized and non-recursive setting in 
[19], [20]. These studies consider a deterministic evolution of the leak intensity and use a continuous 
observation model. We are not aware of any paper that models the time-varying intensity as a random 
process. Most papers study that case when the medium is infinite or semi-infinite since the diffusion 
equation has a closed form solution in that case [15], [17]. The medium model assumed in this paper is 
more general. We also remark that we can extend the results to non-rectangular geometries by using the 
Galerkin approximation [21]. 

While centralized recursive source localization has recieved much interest [15], [17], [20], [22] there 
are very few papers that discuss a distributed solution. A recursive and distributed solution to the problem 
in a Bayesian setting is discussed in [16]. A related paper is [23] that deals with the problem of estimating 
the diffusion coefficient in a distributed and recursive manner. We are not aware of any prior work that 
solves the source localization problem using a distributed and recursive approach in a non-Bayesian 
setting. 

B. Approach 

We show in Appendix |B] that by using Green's technique to solve differential equations it is possible 
to obtain a state-space description for each sensor's observation process. We can then use the IRPE 
algorithm to estimate the iterates in a distributed and recursive manner. 

C. Numerical results 

We use h = h = 100 and diffusion coeffient v = \. The actual location of the source is x* = (37, 48). 
The initial intensity value is taken to be 100, p is fixed at 0.99 and the variance of S{k) is fixed at 10. 
A network of 27 sensors is deployed. To ensure complete coverage of the sensing area, we first placed 9 
sensors on a grid and then randomly deployed 2 sensors in the immediate neighborhood of each of the 
9 sensors. The network is shown in Fig. |2] 
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Fig. 2. A network of 27 sensors. The circles denote the cluster heads and the squares denote the sensors. The source is 
represented by a dot. The arrows indicate the order in which the iterates are passed in the hybrid IRPE. 
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Fig. 3. Estimate of the x-coordinate generated by the standard IRPE, hybrid IRPE and the centralized RPE. Observe that the 
standard IRPE iterates get caught in a local minimum. 



The sampling interval is 10 time units and the measurement noise variance is set to 0.1. In deriving 
the state-space representation, we use ni = n2 = 15. We performed three simulation experiments. 

1) Standard IRPE: 1000 iterations of the IRPE algorithm are used to estimate the source location. 
As discussed, the IRPE algorithm does not use the information that the sensors observe the same 
under-lying process through different observation matrices Hi. 

2) Hybrid IRPE: The network is divided into 9 clusters of size 3. The cluster heads are the sensors on 
the grid. At the beginning of each slot, a cluster head collects the measurements from the sensors 
in the cluster. To estimate the sensor location, 1000 iterations of IRPE are run between the cluster 
heads. In determining the predictor family for each cluster's observations the information that all 
the sensors in the cluster observe the same underlying process is used. But, in the inter-cluster 



19 



60 
55 

"E 

>-45 
I 40 



100 200 300 400 500 600 700 800 900 1000 
Iterations 

Fig. 4. Estimate of the j/-coordinate generated by the standard IRPE, hybrid IRPE and the centralized RPE. Observe that the 
standard IRPE iterates get caught in a local minimum. 

processing through the IRPE algorithm this information is not used. 
3) Centralized RPE: All sensors immediately communicate their measurements to a fusion center. In 

this case the information is completely used. The fusion center runs 1000 iterations. 
The results are plotted in Figs. |3] and ID As expected the centralized RPE performs the best. However, 
what is interesting to note is that the standard IRPE does not converge to the correct solution but is 
caught in a local minimum. We also observed this in other simulation runs. However, when the sensors 
are clustered the iterates converge to the correct location. 

IX. Conclusions 

Linear state-space models arise naturally, or as linear approximations to non-linear state-space models, 
in many applications. In an inference setting where the aim is to estimate a quantity of interest, it is quite 
natural for the state-space models to be parametrized by the unknown quantity of interest. In a control 
system setting where the aim is to control the process, it is common to use such incomplete state-space 
models as 'grey-box' descriptions of the system that is to be controlled. Thus, the problem addressed in 
this paper is important in both of these settings. 

In Section IVIII we could only give a qualitative description of the tradeoff between centralized and 
distributed schemes. It is therefore of interest to find good bounds on the performance of the IRPE 
and RPE schemes that can be used to quantify the loss in performance. Also, to truly understand the 
performance of the algorithm in practical settings, we need to obtain convergence results when there 
are communication errors. Further, we have considered a simple class of networks where the topology 
is fixed. It is important to obtain an algorithm that is similar to the IRPE for networks with a random 
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and time-varying topologies. Finally, as mentioned in Section IVIII the result and analysis extend easily 
to the case where there is an open-loop input to the system. An important extension is to obtain similar 
convergence results for some common classes of closed-loop inputs using the techniques discussed in 
Appendix 7. A of [1]. 
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Appendix 

A. Proof of Theorem |2] 

We take the following approach to prove Theorem |2] First, we consider a centralized system where 
the sensors immediately communicate all their measurements to a fusion center. For this system we use 
the RPE algorithm to generate a sequence of iterates, and then we show that this iterate sequence is 
identical to the iterate sequence generated by the IRPE algorithm. We complete the proof by proving 
that the iterates converge to a local minimum of f{x). 

In what follows, we extensively use the notion of block vectors and matrices. For positive integers a 
and 6, let TWaxb be the vector space of all real matrices of dimensions a x 6. A block vector in Maxh 
a vector whose elements are from M.axh- The length of a block vector is the number of block elements. 
In a similar manner, block matrices in Maxb are matrices where each element is itself a matrix from 
M^axb- While writing block matrices we will allow for a slight abuse of notation and use and / to 
denote the zero and identity matrices, respectively. Their dimensions can be unambiguously fixed from 
the dimensions of the other blocks in the block matrix. We will use U^, b < m, to denote the unit block 
vector in -Max a of length m, with the b-th block equal to the identity matrix in A^^xa and all the other 
blocks equal to the zero matrix in Maxa- 

We allow i,j to take values in the set X = {1, ... , m}. We define 6[-] as the Kronecker delta. Recall 
that the dimension of the matrices Qi{k;x) is q, the dimension of the measurement rj(fc) is p, and the 
dimension of the parameter vector x is d. Also, recall that for any random process {Y{k;x)} that is 
parametrized by x, y{k) denotes the sample path of Y{k\x*). 

1 ) State-space model for sensor observations: Without loss of generality, assume that each time slot 
has duration of m time units. Consider a hypothetical centralized scheme where at time mk + j, sensor j 
communicates rj{k+l) to the fusion center over a perfect delayless link. For i / j, sensor i communicates 
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a predetermined constant value, say 0, that does not convey any information about the value taken by 
the parameter x. 

Denote the sequence communicated by a sensor i by {fi{mk + j)}, with 



ri{mk + j) = ri{k + l)5[i - j]. 
Next, denote the observation sequence at the fusion center by {r{mk + j)}, where 



(25) 



(26) 



r{mk + j) = [ri{mk + j) . . . Vmimk + j)] =Ujrj{k + l). 

We now consider the problem of estimating x from observation sequence {r{mk + j)}. We show that 
the random process R{mk + j;x) can be represented as the output vector of a suitably defined state- 
space system. For this, we first use the relations in ([Hi to obtain the equations describing the evolution 
{Ri{mk + j;x)}. Note that from (|25] ). we have 



Ri{mk + j; x) = Ri{k + l;x) 6[i - j]. 
Let Di{x) be the following ui X fn block matrix in .Adqxq' 



Di(x) 



/ 

/ 



D,{x) 






/ 





Observe that 



when j / 1 
Ul^Diix) whenj = l. 



Also, define Hi = Hi (C/f )^ , and note that 



Him = H, {u\fm = H,5[j-i]. 



Define 9^(0; x) = Ul 9,(0; x), and 



@i{mk + j;x) 



Utj_^,ei{k + l;x) if j<i 



^ rn+l-(j-i) 



Qi{k + 2; x) if j > i, 
W^imk+j;x) = Ul W^{k + l;x) d[i - j], 
Vi{mk+j) = V^{k+l) d[i-j]. 



(27) 



(28) 



(29) 



(30) 



(31) 



(32) 
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The following is an illustration for Qi{mk + j; x) with i = 3 and j = 2, 3, 4: 

mA; + 2 mA; + 3 mk + 4 





e3(A; + l;x) 







e3(A; + l;x) 














Q3ik + 2-x) 



Claim 1: For all n > 0, we have 

Qi{n + l;x) = Di{x)@i{n; x) + Wi{n; x), 
Ri{n + l]x) = HiQiin+l^x) + Vi{n + l). 



(33) 
(34) 



Proof: Let n = mk + j. Substitute for Q)i{n;x) and from (|3TI ) in the right hand side 

(RHS) of (HD. For i < j, from ^ we obtain 



RHS of (131 = Di{x)Ul_.^^ @i{k + 1; x) + 
= Ul.ei{k + l;x) 



For i = j, using 



^Q^{k + l■,x) 
= Qi{mk + j + 

and ([Hi, we obtain 

RHS of dH = Di{x)U\ Q^{k + 1; x) + ?7^, Wi(A; + 1; x) 
= Ulfii{x)Qi{k + l-x) + W,{k + 1; x) 
= Ul{D,{x)e,{k + + Wi{k + 1; x)) 
= J7^ei(fc + 2;x) 
= Qi{mk + j + l;x). 



24 



Finally, when i < j, from ( 1291 ) we have 

RHS of (131 = Di{x)Ul_^._.^^^ G,{k + 2;x) + 

= f/L-O+i-.)-i0^(^ + 2;^) 
= Qi{mk + j + l;x), 



thus, completing the proof of the relation in (1331 ). 

We next prove the relation in (l34l ). At first, we consider the case when j ^ m and show that the 
following relation holds: 

HSiimk + j + l;x) + Vi{mk + j + l) = R,{k + l)6[i - j - 1]. 

Let i 7^ J ■ + 1, and note that from (l32l ) we have Vi(?^ + 1) = 0. Furthermore, from the definition of 
Qi{n;x) in (13T1 ) we obtain 



i?ie,(mfc + (j + l);2:) = 



[/t,e,(fc + l;x) ^>J + l 



Using the expression in (l30l ). we can immediately verify that HiQi{mk + (j + =0. Therefore, 

we have Ri{n + 1; x) = for i 7^ j + 1. 

When i = j + 1, Vi(?nA: + j + 1) = yj+i(A; + 1) and 

Hj+iQi{mk + j + l-x)= Hj+iUl = Hj+iQj+i{k + 1; x). 

Therefore, from (HJ we see that Ri{mk + j + \;x) = Ri{k + l; x) for i = j ' + 1, thus, concluding the 
proof of relation (l34b for j in. 

When j = m, by using arguments similar to that of the preceding case j / m, we can show that 

H.Qiimk + j + 1; x) + Vi{mk + j + I) = Ri{k + 2)6[i - 1], 

thus, completing the proof. ■ 
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Now, by combining the equations in (l33])-(l34l) for i G T, we provide evolution equations for {R{n; x)}. 
Define 

F{x) =diag{Fi{x),...,Fm{x)) , 
H{x)=dmg {Hi{x),...,H,^{x)), 
Qi{n;x) 



6(n; x) 



W{n; x) 



@m{n;x) 

Using the relations in (1331) and (134] ). we can write 



Wi{n; x) 
Wmin;x) 



V{n; x) 



Vi (re; x) 
Vm{n;x) 



(35) 



e(n + l;x) = D{x)e{n; x) + W{n; x), 
R{n + l;x) = HQiin + 1; x) + V{n + 1). 



(36) 
(37) 



To apply the RPE we need to determine a predictor family for ^(re; x*) that is parametrized by x and 
is asymptotically optimal at x = x*. We do this in the next section. 

2 ) Time-Invariant Kalman Predictor for Centralized System: Let us first obtain the time-invariant 
Kaknan predictor for Ri{n] x). Fix re = mk + j and define 

(^i,fe+i(x;rf) ifi<i 



(38) 



9i,n{x'^'^i ) = 9i,k+i{x;ri) 6[j - i]. 



(39) 



Note that 4>i,k+i{x*; rf) and gi^k+i{x*',r^), are the time-invariant Kalman predictors for 0j(/c + l; x*) and 
Ri{k-\-l; X*), respectively. Therefore, from (|3T] ) we can conclude that 4'i,nix*; f"~^)), resp. gi,n{x*;f^~^), 
is asymptotically optimal for Qi{n;x*), resp. Ri{n;x*). 
Define 



Gi{x) = Ul, G,ix) 



(40) 



and Fi{x) = Di{x) — Gi{x)Hi. The matrix Fi{x) will have the same form as Di[x) in (1281) but with 
Di{x) replaced by Fi{x). Similar to Claim [H we can show that 



(/)i,„+i(x;r") = Fi{x)4ii^n{x;r^ ^) + Gj(x)rj(re), 
5i,„+i(x;f") = Hi{x)(i)i^n{x;fi). 



(41) 
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We can immediately obtain a predictor family for @{n;x*) and {Rn{x*)} that is asymptotically optimal 
at X = X* as follows: 



-n-l^ 
m > 



= [9i,n(2;;rJ'-i)...5„,„(j;;f;^-i)]^. 
Furthermore, from (|4TI ) one can verify that 

<^„+i(rE;f") = + G(rE)f(n), 

5„+i(a;;f") =#<^„+i(x;f'^), (42) 

where 

= diag . . . , . (43) 

3) RPE Algorithm for Centralized System: Here, we use the RPE algorithm to estimate x from {f(n)}. 
As we mentioned earlier, r{n) contans the same information as f(n) about the true value of x. Define 
for £ = 1 , . . . , d, 



vWf( 



dF{x) 

dxW 



v(^)G(x) 



dGjx) 

dxW 



We define the iterates {x„} as follows: 



hn+l 




H 







1pn+l 









H 




An+1 



(^) 
■n+1 



f(n + 1) - /i„. 



+1) 



X, 



n 



+ 

r 



i£n+l • • • 2in+l 







An+2 





F{Xn+l) 

vWf(x„+i) F{x n+1) 



+ 



G{Xn+l) 

VWG(x„+i) 



(44) 

(45) 
(46) 
(47) 
(48) 

f(n + l). (49) 



Here, a{n) = a^+i for n = nik + j for j = 1, . . . ,m — 1. Next, we assign the initial values for the 
recursion. Recall that the IRPE algorithm in ([T< 
all i and i, and xq = Xs- We let xq = Xs, and 



recursion. Recall that the IRPE algorithm in (fT9l) is initialized with the values tpi^i = tpi^s, £,fi = ifl 











m _ 

' so — 















(50) 
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where 'ipi^s = and = t/^.^-'^j for all i and I. 

4) Rest of the proof: Finally, here we show that Xn = Zj,k+i- Recall that ^j^^ and xTk generated 
in the IRPE algorithm ([14ll-([T9ll. Define for / = 1, . . . , d, 



i+1 '^i,k+l 



if j <i 



(51) 



(52) 



if j < i 

We next establish the following lemma. Once we prove this lemma we can use induction to show that 
the iterates generated by the IRPE algorithm are the same as those generated by the centralized scheme. 



Lemma 1: Let n = mk + j. If Xn = Zj^j^i and 



An+1 



Al,n+1 ■ • • /'^m,n+l 



T 



for £ = !,..., d, 



(53) 
(54) 



then Xn = Zj+i^fe+i, and 



+2 



X.n+2 



[^^1,71+2 ■ • • 1pm,n+2] ' 

A{e) M) 

Al,n+2 • ■ ■ A7Ti,n+2 



for t=l,...,d. 

Proof: For a block vector A, let A^'^'^ denote its i-th block. Substituting for ^mk+j+i from (1531 ) in 
(l44l ). and noting from (|35] ) that H is a. block diagonal matrix with the (i, f)-th block equal to Hi, we can 
see that 

Using the definition of -ipi^rnk+j+i from (ISTl ) and noting from (l30l ) that HiU'j = Hi5[j — 1], we obtain 

if i > j + 1 

HiUlipi^k+i ifi= j + l 

if i < j + 1 
if i ^ j + 1 

iJj+iV'j+i.fc+i if i = j + 1. 
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Using ([T4l ) we replace ifj+i'^j+i.fc+i by hj^i^k+i and write 

Therefore, hmk+j+i = Similarly, we can see that for all £, 



''mfc+j+i ~ i+1 ^j+i,k+r y-'-^' 



Substituting in (1451) for hmk+j+i from above and for fmk+j+i from (1261) we get 

Emfc+i+l = - rj+i{k + 1)). 

Observe that ej+i^fc+i = /ij+i^fc+i — rj^i{k + 1) from (fTSl) . so that 

^mk+j+l 

Since Xmk+j = -^i./fc+i it follows that x^J^i^_^_j = zf\+i- Substituting from (1551 ) and (l56l ) in (l46l) we get 
for all I, 

_ (^) 

The last step follows from ([T6l ). Therefore, from (l47l) and (fTTl) we can conclude that x^^j^^j^i = z^j^i 
and from ( |48] ) and ( fTSl ) that Xm^+j+i = z^+i fc+i. This completes the first part of the proof. 

Let us next consider the case when j G {1, . . . , m — 1}. In ( |49l ) let us replace Xn+i with Zj+i fc+i. 
Note from (l35l) . respectively (1431 ). that ^(zj+i fc+i), respectively (^(zj+i^fc+i), is a block diagonal matrix 
with (i,i)-th block equal to fc+i), respectively Gi{zjj^i^f^j^i). Substituting for V'n+i from ([53]) 

and f(n + 1) from (l26l) in ( |49l) we can write 

=-^i(2j+i,fc+i)'0j,mfc+i+i + Gi{zj+i^k+i)ri{mk + j + 1). 

Let us substitute for G'j(zj_|_i^fc-|_i) from ( |43] ). for ■0j „^fc_|_j_|_i from dST]) and for fi{mk + j + 1) from 
Using ( [291 ) we get for i > j + 1, 

V5j^i+j+2 = Fi{zi,k+i)Ul__j^i,k+i + 

= u1_j_-^iJi^k+i 

= '4'i,mk+j+2- 



29 



When i = j + 1, from ([T9l ) we obtain 

Fi{zi,k+i)tpi,k+i + U'^mGi{zi^k+i)riik + 1) 

= i^i,mk+j+2- 



When z < j + 1, we have 



= Fi{z,^k+l)Ul_^._.^'ll^i,k+2 + 



= V'i,mfc+j+2- 

Equation (|53]) can be shown similarly for the case j = m. The proof of relation (l54l ) is very similar to 
that of relation (1531 ) and therefore, it is omitted. ■ 
Observe that the initial values for the RPE algorithm in ( [49l ) and the IRPE algorithm in ( fT9l ) are 
choosen such that xq = zo,i; and 

■01 = [^1,1 • • • i'm.lf , 

^ for all 



Xi 



Al,l • • • Xm,l 



By using Lemma [U and the induction on k, we can conclude that Xmk+i = Zi^k+i for all A; > 1 and 
« G X. To complete the proof, we only need to show that the sequence converges to a minimum of 
f{x), which is done in the following. 

Lemma 2: The sequence {x„ } generated by (i44\l—l[49\) converges to a local minimum of the function 
f{x), defined in ( |7?1 ). over the set X w.p.l. 

Proof: By the assumptions of Theorem |2j it follows that the conditions of Theorem [T] are satisfied. 
Therefore, by Theorem \T\ the iterates x„ converge to a local minimum over the set X of the following 
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function 



1 ^ 

f(x) = lim — > E 



TV 



TV 



n=l 



n=l j=l 



R{n;x*)-gnix,R''ix*)) 

n=l j=l 
^ N m 

I — ^■oo iV — ^ — 7 L 



B. State-space Model for the Sensor Measurements 

A standard technique to solve partial differential equations with boundary and initial conditions is to 
use Green's function. For the equation in (l24l ) the solution can be written as follows 



C{y,t;x) 



I{t)5{z — x)Gf{y, z,t — T)dz dr. 



'0 JD 

Here, r and ^ = (^1,22) are parameters of integration, and Gj is the following Green's function 

2 



1 2 

ni,n2=l 1=1 



li 



COS ; cos 



Evaluating C{y,t;x), we obtain 



^ ^ ni=ln2=li=l ^ * 



C{y,t]x) 



where 



niT^yi \ ( riiTTXi 
cos 



(57) 



Pn„n2 = exp ^-l/TT^ ^ . 

To get a convenient approximation we will use a sufficiently large, but fixed number of terms in the 
convergent infinite series in (l57l) . We will let rij, i = 1,2, vary from 1 to fii, where the integers > 
are chosen large enough to provide a sufficiently good approximation. Therefore, 

ft "1 "2 2 ^2 /^..^„,.\ / rt 

k _ 

(58) 



C{y,t;x] 



I ri ni na 2 / 2 

^ ^ ni=ln2=li=l 
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Define 



0oW = ^^*At) dr. 

2 

P{y,ni,n2) = JJcos 



2 



k 



A'{x, ni, n2) = Y\ cos ^ ^'^^' ^ , 



With this notation we can write (1581 ) as 

fii n2 

C{y,t;x) = Q'oit) + E 0'n„n.(i; '^i, "2). (59) 

ni=l 712=1 

From the function /(t) in (l22l ). we have 

e'o(A: + l) =e[)(A;)+/(A: + l), 
e'„,^,„^(A; + 1; =f3n„nA,n.{k; x) + m, n2)I(A; + 1) f ^ " ^ " 
Furthermore, define 

7ni,n2 = 1 + n2(ni - 1) + n2. 
From ( [58] ) and ( [59l ) we can write C{y, k + l;x) as the output of the following state-space system: 

e'{k + 1; x) =D'Q'{k; x) + B'{x)I{k + 1) 
C{y,k + l;x)=H'{y)e'(k + l;x), (60) 

where 

©Wi =[^o(A: + 1)^1,1 (A: + 1) . . . ^i.^^CA: + 1) 02,i(A: + 1) . . . ^^^^^^ (A: + 1)]^ , 
is the diagonal matrix with 

Z)'(l,l) = 1, D\^n,,n,nn,,n,)=fiZ,n,, 

B'{x) is the column vector 

= 1 

A'(x,ni,n2)/?!r « -1 
(7„„„2, 1) = Vlr/^ T ' 
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and H'{y) is the row vector with 

H'iy) = [1 P{y, 1, 1) Piy, 1,2)... P(y, 1, na) P(y, 2, 1) . . . P{y, m,n2)] . 
From dH]), (inil and ([2311 we have 

1 

e'(A;+ 
I{k + 1) 

Thus, the system dynamics are modeled using a state-space model similar to ([U that is parametrized by 
the unknown source location x. Notice that in this case 

^ e'{k;x) 



e'{k + l;x) 




D' pB\x) 




@'{k;x) 


+ 


" B'{x) ' 


I{k + 1) 




. P . 




m 




1 



Ri{k + l;x)= H'{si) 



+ Ni(k + 1). 



Qi{k; x) = &j{k; x) 



m 



Hence, complete information is available about the joint statistics of Qi{k]x) and @j{k;x). If the sensors 
do not sense synchronously at the beginning of the slot, the model is more complex. Nevertheless, we 
can still identify a state vector and write the sensor measurements as the output process of a state-space 
system. 



